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One-dimensional  molecular  dynamics  simulations  of  the  onset  of  detonation  have  been  performed 
using  three-body  potentials  which  accurately  reproduce  the  effects  of  endothermic  bond  breaking 
and  exothermic  bond  recombinations.  A  stable  del-nation  wave  of  reasonable  velocity  can  be 
produced  from  the  impact  of  a  plate  upon  a  stationary  array  of  diatomic  molecules.  The  onset  of 
detonation  has  been  studied  using  a  variety  of  potential  forms  including  LEPS  and  Tersoff 
potentials.  Properties  of  the  detonation  front,  including  front  velocity,  reaction  zone  width, 
ar,d  product  distribution,  have  been  studied  as  functions  of  barrier  height,  oxothermicity,  and 
other  relevant  potential  function  parameters.  A  model  is  proposed  to  explain  the  sel f -regulating 
detonation  front  velocities  observed  in  these  simulations  in  terms  of  atomic-scale  kinetics. 

].  INTRODUCTION  functions  which  are  capable  of  accurately 

The  study  of  the  onset  of  detonation  at  the  describing  bond-breaking  and  bond-forming  behavior 

molecular  level  is  a  difficult  experimental  task  in  simple  model  systems  of  detonating  solids.  We 

because  of  the  short  time  scales  involved  and  the  demonstrate  that  one-dimensional  molecular 

destructive  nature  of  the  detonation  process.  dynamics  simulations  using  these  potentials 

Computer  simulations  are  an  obvious  alternative,  produce  self-sustaining  detonation  waves  with 

but  the  continuum  methods  which  are  so  valuable  stable  velocities  independent  of  initiation 

for  predictions  at  the  macroscopic  level  are  conditions.  Finally,  we  offer  a  molecular 

unable  to  give  information  regarding  the  primary  interpretation  for  the  stabil ity  of  the  detonation 

chemical  processes  occurring  in  the  reacting  velocities, 

solid.  Molecular  dynamics  simulations,  on  the 

other  hand,  are  feasible  under  precisely  the  2.  THE  LEPS  POTENTIAL 

conditions  of  interest:  short  times  and  relatively  One  form  of  potential  function  which  can  give 

small  numbers  of  molecules.  Indeed,  molecular  qualitatively  correct  bond-breaking  and  bond- 

dynamics  simulations  of  nonreactive  shock  waves  forming  behavior  for  simple  systems  is  the  LEPS 

have  been  performed  for  more  than  two  decades.  (London,  Eyring,  Polanyi,  Sato)  potential, 

For  the  study  of  detonation,  however,  the  developed  in  the  1930' s  to  describe  the  H3 

simulations  must  explicitly  include  the  chemical  potential  energy  surface.  In  the  form  devised  by 

reactions  which  provide  the  energy  to  sustain  the  Tully1,  a  three-body  LEPS  potential  function  can 

process.  This  requires  the  use  of  many-body  be  described  in  terms  of  2x2  matrices  as  follows: 

potential  functions  which  correctly  reproduce  the  For  each  pair  of  atoms  ij  in  the  three-body 

sal ient  features  of  detonation  reactions,  but  yet  system,  a  diagonal  2x2  matrix  is  constructed  in 

are  simple  enough  to  make  the  simulation  which  Hu  -  Vbond(ru)  and  H22  =  Vnonbond(rtJ),  where 

computationally  fca'ible.  ru  is  the  scalar  distance  between  atoms  i  and  j, 

We  report  two  types  of  many-body  potential  is  a  typical  diatomic  bonding  potential  such 
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The  LEPS  potential  is  particularly  suitable 
for  one-dimensional  simulations  because  in  such 
systems  each  atom  has  exactly  two  nearest 
neighbors  and  therefore  forms  the  central  atom  of 
a  three-body  system.  The  total  potential  energy 
of  a  one-dimensional  crystal  can  be  expressed  as 
a  sum  of  overlapping  LEPS  three-body  potentials 
in  which  each  atom  (except  those  on  the  ends) 
appears  as  a  left-hand,  central,  and  right-hand 
member  of  three  different  trimers.  The  potential 
energy  and  atomic  forces  can  be  expressed 
analytically,  so  the  method  is  computationally 
efficient. 

We  have  recently  reported2  preliminary  results 
of  some  one-dimensional  molecular  dynamics 
simulations  using  LEPS  potentials  designed  to 
represent  the  features  of  detonating  nitric  oxide 
system,  2N0  -  112  +  02.  The  initial  NO  crystal  is 
arranged  in  the  form  -  .  NO  ON  NO  ON  • • •  so  that 
formation  of  the  appropriate  products  is  possible. 
The  NNO  and  1100  LEPS  potentials  are  parameterized 
to  match  experimental  exothermicities  and 
calculated  barrier  heights  for  the  three-body 
reactions  N  +  110  -  N2  +  0  and  N  +  02  -•  NO  +  0.  The 
detonation  is  initiated  by  impact  of  a  small 
nitric  oxide  plate  upon  one  end  of  the  larger 
crystal.  The  simulation  can  readily  be  followed 
over  the  course  of  20  picoseconds,  during  which 
time  the  detonation  front  propagates  through 
several  thousand  atoms. 

The  most  striking  feature  of  these  simulations 
is  that  the  detonation  velocity  is  Independent  of 
the  plate  impact  velocity,  as  long  as  the  impact 


injfi  '<■  >  cj  ■  :  .  :  :  ruptures  which 
initiate  the-  c:;.i  I  . •  ..  ;  : ,  -,s  .  i  i  the  plate 

initially  overdriven,  but  (...Kyiv  returns  to  Us 
steady-state  value.  !f  the  impact  velocity  is  too 
low  to  initiate  detonation,  it  a  gradually 
decaying  ncnreact  ivo  •  r<  •.  *  w  y . <  ‘  u  I  rop-cat- 
through  the  system.  The  iacx  ot  dependence  o'  the 
detonation  velocity  on  initiation  conditions  is 
in  agreement  with  experimental  observations  of 
condensed- phase  detonation. 

figure  1  shows  the  atomic  density  relative  to 
the  front  for  a  plate  impact  speed  of  6  km/sec, 
averaged  over  the  course  of  a  20  picosecond 
simulation.  Most  of  the  chemical  reactions  take 
place  in  the  high-density  region  within  400 
Angstroms  of  the  front;  behind  this  region  the 
density  falls  off  rapidly  as  the  products  and 
unreacted  fragments  separate.  Further  reaction 
in  this  gas-phase-like  region  is  less  likely 
because  of  the  larger  average  interatomic  spacing 
and  small  relative  velocities  which  prevail  in 
this  region. 


Position  (A)  relative  to  front 


FIGURE  1 

Density  (g/mqle/A)  versus  position  relative  to  the 
front  for  a  plate  impact  speed  of  6  km/sec. 
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One  difficulty  with  the  ;Er$  (iotr-nti.il  ■■  tn.> 
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the  method  to  higher  dimensions  with  arbitrary 
numbers  of  near  neighbors.  We  have  therefore 
investigated  the  use  of  another  form  of  empirical 
reactive  potential  based  on  tre  work  of  In 
The  Tersoff  method  uses  a  modified  pair  potential 
in  which  the  bonding  part  of  the  interaction 
between  a  pair  of  atoms  i  and  j  is  weakened  by  an 
amount  proportional  to  the  number  of  other  near 
neighbors  which  compete  for  bonding  to  atom  i, 
i.e., 

E  -  l  E, 
where  1 

E,  -  d/2)I  (V„(rtJ)  -  B1J-VA(rlJ)) 

Here  V-^  represents  a  repulsive  potential,  V.  is  an 
attractive  potential,  and  B,_,  represents  a  many- 
body  coupling  between  bond  ij  the  local 
coordination  of  atom  i.  These  three  functions  can 
be  implemented  in  a  variety  of  ways‘-i,  giving  the 
Tersoff  potential  a  great  deal  of  flexibility. 
In  particular,  the  coupling  term  Bu  may  depend  on 
the  angle  of  a  competing  bond  ik  to  the  bond  ij, 
thereby  incorporating  directed  valence  effects  of 
importance  for  covalent  compounds  in  two  and  three 
dimensions. 

We  have  constructed  a  one-dimensional  Tersoff 
potential  to  represent  the  generic  reaction  2AB 
-  Aj  +  Bj,  where  the  AB  molecules  have  been  given 
a  bond  strength  of  4  eV  and  the  A2  and  B2 
molecules  have  each  been  assigned  a  bond  energy 
of  5  eV  so  that  the  overall  reaction  is 
exothermic.  An  initial  AB  molecular  crystal  was 
constructed  in  a  manner  identical  to  the  nitric 
oxide  simulations  described  above,  •••  AB  BA  AB 
•••,  and  the  stability  of  the  reactant  solid  was 
ensured  by  incorporating  barrier  heights  of  0.5 
eV  for  the  two  possible  three-body  reactions  A  + 
AB  -  A2  +  B  and  B  +  BA  -  B2  +  A  (which  are 
energetically  equivalent  in  this  case)  and  weak 
long-range  van  der  Waals  attractions  between  AB 
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described  above,  molecular  dynamics  simulations 
of  the  detonation  of  this  system  were  carried  out. 
with  the  detonation  initiated  by  •.-pact  .it  a  -.ass 
of  variant.;-  velocity  cn  tr-e  reactant 
Although  the  density,  exothermicitv,  and  the 
potential  surfaces  in  this  case  were  quite 
different  from  those  employed  in  the  nitric  oxide 
simulations,  the  results  were  qualitatively 
similar.  For  impact  velocities  ranging  from  4 
km/ser  to  9  km/sec,  the  detonation  velocity 
rapidly  stabilized  to  a  constant  value  of  12.6 
km/sec.  Exothermic  chemical  reactions  in  a  high- 
density  region  behind  the  front  provided  the 
energy  to  sustain  the  detonation.  Impact 
velocities  less  than  4  km/sec  or  greater  than  9 
km/sec  did  not  produce  detonation;  in  these  cases 
the  disturbance  propagated  through  the  chain  as 
a  gradually  dissipating  nonreactive  shock  wave. 

4.  VELOCITY  REGULATION  MECHANISM 

The  self-limiting  detonation  velocity  observed 
in  the  simulations  can  be  understood  by 
investigating  the  behavior  of  an  isolated  three- 
body  reaction  process,  A  +  AB  -  A2  +  B.  Classical 
trajectory  calculations  on  the  Tersoff  potential 
surface  described  above  lead  to  the  results  shown 
in  Figure  i,  where  the  velocity  of  product  atom 
B  is  plotted  as  a  function  of  the  velocity  given 
to  an  incident  reactant  atom  A.  The  straight  line 
marks  the  condition  that  the  velocity  of  the 
outgoing  atom  equals  the  velocity  of  the  incoming 
atom.  The  lower  value  of  the  curve  is  determined 
by  the  smallest  velocity  of  atom.  A.  which  ran 
overcome  the  potential  barrier  to  reaction.  The 
upper  value  of  the  curve  represents  the  point 
beyond  which  the  impact  of  atom  A  leads  only  to 
endothermic  fragmentation  (A  +  A  +  B)  or  to  no 
reaction  (A  +  AB). 

If  detonation  in  one  dimension  is  viewed  as 
being  driven  by  a  series  of  three-body  reactions 
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FIGURE  2 

Outgoing  velocity  of  atom  B  following  the  reaction 
A  +  AB  -  A,  4  B,  as  a  function  of  the  velocity  of 
the  incoming  atom  A  relative  to  a  stationary  A6 
molecule.  The  straight  line  indicates  equal 
incoming  and  outgoing  velocities. 

in  which  the  products  of  one  reaction  col  1  ide  with 
a  molecule  to  initiate  a  new  reaction,  then  the 
detonation  regime  will  be  limited  to  the  range  of 
impact  velocities  which  cause  initiation  of  the 
exothermic  reactions  A  +  AB  -  Aj  +  B  and  B  4  BA  - 
Bj  4  A.  This  analysis  should  also  be  valid  in 
higher  dimensions,  as  cross  sections  for  bond 
dissociation  generally  increase  for  increasing 
atomic  collision  energies. 

For  incoming  velocities  near  the  low  end  of 
the  reaction  range,  the  velocity  of  the  product 
atom  is  greater  than  that  of  the  reactant  atom. 
This  is  due  to  the  net  exotherraicity  of  the 
reaction,  and  has  the  effect  of  accelerating  the 
atomic  velocities.  For  higher  incident 
velocities,  however,  the  outgoing  atom  has  a  lower 
velocity  than  the  incoming  atom.  This  occurs 
because  the  trajectory  interacts  with  the 
repulsive  wall  of  the  potential  surface, 
increasing  the  amount  of  energy  deposited  in  the 
internal  vibrational  mode  of  the  product  molecule 
Aj.  If  a  series  of  isolated  reactions  in  the 
orientation  present  in  the  chain  described  above 
are  initiated  with  an  atom  anywhere  within  the 
range  of  reacting  velocities,  the  series  will 
"self  regulate"  to  the  point  at  which  the  two 
curves  cross  in  Figure  2.  Of  course,  in  the  solid 


the  detonation  velocity  will 

iigm :  leant  Iy  higher  than  the  above  anal., 

;  .  because  the  detonation  wave 

collective  disturbance  which  can  propagate  faster 
than  individual  atomic  velocities.  This  effect 
can  tv-  accounted  for  quantitatively  by  incluqir-: 
u.  t  •  analysis  the  effect  of  the  density  •; : 
chain,  or  more  accurately  the  fraction  of  the  i.;.n 
cell  length  excluded  by  repulsive  atomic  core 
potential s . 

5.  CONCLUSION 

Molecular  dynamics  simulations  of  detonation 
can  be  carried  out  in  one  dimension  with  many-body 
potential  functions  which  incorporate  realistic 
reaction  dynamics.  These  simulations  can  provide 
valuable  insight  into  the  chemical  mechanisms 
underlying  the  initiation  and  propagation  of 
detonation.  The  Tersoff  potential  is  particularly 
well  suited  to  the  extension  of  these 
investigations  Vo  two  and  three  dimensions. 
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